Открытый курс по машинному обучению. Сессия № 3

Автор материала: Дмитрий Сергеев, Zeptolab. Материал распространяется на условиях лицензии Creative Commons CC BY-NC-SA 4.0. Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала.

Домашнее задание № 9

Временные ряды

Заполните код в клетках и выберите ответы в веб-форме.

Импортируем необходимые библиотеки


In [83]:
import warnings
warnings.filterwarnings('ignore')

import random
random.seed(42)

import seaborn as sns
import pandas as pd
import numpy as np

from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import TimeSeriesSplit, cross_val_score

import matplotlib.pyplot as plt
%matplotlib inline

def mean_absolute_percentage_error(y_true, y_pred): 
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

Для работы возьмем реальные часовые данные по суммарному просмотру рекламы в одном из приложений.


In [84]:
df = pd.read_csv('../../mlcourse_open//data/ads_hour.csv',index_col=['Date'], 
                 parse_dates=['Date'])

In [85]:
with plt.style.context('bmh'):    
    plt.figure(figsize=(15, 8))
    plt.title('Ads watched (hour ticks)')
    plt.plot(df.ads);


Тут есть и тренд, и различная сезонность, и немножко выбросов.

В этом домашнем задании мы не будем пытаться подобрать наилучший ARIMA-процесс, порождающий наш ряд, и не затронем библиотеку Prophet от Facebook. Работа с этими методами очень подробна расписана во второй и третьей части ноутбуков по временным рядам в репозитории открытого курса по машинному обучению.

А здесь мы сосредоточимся на построении и отборе признаков для временных рядов.

Возьмем две полезные функции из лекции – для расчета среднего значения по таргету и для подготовки обучающего и тестового датасетов:


In [86]:
def code_mean(data, cat_feature, real_feature):
    """
    Возвращает словарь, где ключами являются уникальные категории признака cat_feature, 
    а значениями - средние по real_feature
    """
    return dict(data.groupby(cat_feature)[real_feature].mean())

def prepareData(data, lag_start=5, lag_end=14, test_size=0.15, target_encoding=False):
    """
        series: pd.DataFrame
            dataframe with timeseries

        lag_start: int
            initial step back in time to slice target variable 
            example - lag_start = 1 means that the model 
                      will see yesterday's values to predict today

        lag_end: int
            final step back in time to slice target variable
            example - lag_end = 4 means that the model 
                      will see up to 4 days back in time to predict today

        test_size: float
            size of the test dataset after train/test split as percentage of dataset

        target_encoding: boolean
            if True - add target averages to the dataset
        
    """
    data = pd.DataFrame(data.copy())
    data.columns = ["y"]
    
    # считаем индекс в датафрейме, после которого начинается тестовыый отрезок
    test_index = int(len(data) * (1 - test_size))
    
    # добавляем лаги исходного ряда в качестве признаков
    for i in range(lag_start, lag_end):
        data["lag_{}".format(i)] = data.y.shift(i)
        
    
    data.index = data.index.to_datetime()
    data["hour"] = data.index.hour
    data["weekday"] = data.index.weekday
    data['is_weekend'] = data.weekday.isin([5,6])*1
    
    if target_encoding:
        # считаем средние только по тренировочной части, чтобы избежать лика
        data['weekday_average'] = list(map(code_mean(data[:test_index], 'weekday', "y").get, data.weekday))
        data["hour_average"] = list(map(code_mean(data[:test_index], 'hour', "y").get, data.hour))

        # выкидываем закодированные средними признаки 
        data.drop(["hour", "weekday"], axis=1, inplace=True)

    data = data.dropna()
    data = data.reset_index(drop=True)
    
    
    # разбиваем весь датасет на тренировочную и тестовую выборку
    X_train = data.loc[:test_index].drop(["y"], axis=1)
    y_train = data.loc[:test_index]["y"]
    X_test = data.loc[test_index:].drop(["y"], axis=1)
    y_test = data.loc[test_index:]["y"]
    
    return X_train, X_test, y_train, y_test

Воспользуйтесь предложенными функциями и приготовьте из исходного датасета df необходимые для обучения датафреймы. Для тестирования отложите 30% данных, в качестве начального лага возьмите значения временного ряда 12 часов назад, а в качестве последнего лага - значения ряда 48 часов назад. Таким образом модель будет способна строить предсказания на 12 часов вперёд, имея фактические наблюдения за предыдущие полтора дня.

Также отмасштабируйте признаки, для этого воспользуйстесь StandardScaler и создайте переменные X_train_scaled и X_test_scaled. Не забудьте, что обучать scaler нужно только на тренировочной выборке, чтобы информация о средних отклонениях не просочилась из будущего в прошлое.


In [87]:
scaler = StandardScaler()
X_train, X_test, y_train, y_test = prepareData(df, test_size=0.3, lag_start=12, lag_end=48)

In [88]:
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)

Теперь обучите простую линейную регрессию на полученных данных:


In [89]:
lin_reg = LinearRegression()

In [90]:
lin_reg.fit(X_train, y_train)


Out[90]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

Оценим качество работы линейной регрессии на тренировочной части при помощи кросс-валидации. Для этого необходимо сначала создать объект-генератор разбиений для временного ряда при помощи TimeSeriesSplit. Число фолдов задайте равным 5. Затем воспользуйтесь методом cross_val_score, передав туда в качестве параметра cv созданный генератор разбиений. Метрикой качества будет выступать neg_mean_absolute_error.

Не забудьте взять среднее от полученного значения и домножить его на -1.


In [91]:
tscv = TimeSeriesSplit(n_splits=5)
scores = cross_val_score(estimator=lin_reg, cv=tscv, X=X_train, y=y_train, scoring='neg_mean_absolute_error')
-np.average(scores)


Out[91]:
5076.311824059033

Вопрос 1. Чему равно среднее значение MAE на кросс-валидации?

  • 4876
  • 50962872
  • 5076
  • 0.638

Теперь посмотрим на результаты работы линейной регрессии на отложенной выборке. В этом нам поможет функция plotModelResults. В ней Вам нужно дописать часть, ответственную за расчет доверительных интервалов. Для расчета отклонения воспользуйтесь только что реализованным способом оценки качества на кросс-валидации – оттуда понадобится среднее значение MAE и стандартное отклонение этой величины. Доверительные интервалы строятся по формуле $\hat{y} \pm (\text{error} + scale * \text{standard deviation})$


In [92]:
def plotModelResults(model, X_train, X_test, plot_intervals=False):
    """
    Строит график прогнозных и фактических значений, а также доверительных интервалов прогноза
    
    """
    
    # получаем предсказания по модели
    prediction = model.predict(X_test)
    
    plt.figure(figsize=(20, 7))
    plt.plot(prediction, "g", label="prediction", linewidth=2.0)
    plt.plot(y_test.values, label="actual", linewidth=2.0)
    
    if plot_intervals:
        cv = TimeSeriesSplit(n_splits=5)

        scores = cross_val_score(estimator=model, 
                                 cv=tscv, 
                                 X=X_train, 
                                 y=y_train, 
                                 scoring='neg_mean_absolute_error')
        mae = -np.mean(scores)
        deviation = np.std(scores)
        
        
        scale = 1.96
        lower = prediction - (mae + scale*deviation)
        upper = prediction + (mae + scale*deviation)
        
        plt.plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)
        plt.plot(upper, "r--", alpha=0.5)
        pass
        
    mae  = mean_absolute_error(prediction, y_test)
    mape = mean_absolute_percentage_error(prediction, y_test)
    plt.title("MAE {}, MAPE {}%".format(round(mae), round(mape, 2)))
    plt.legend(loc="best")
    plt.grid(True);

Для визуализации коэффициентов модели воспользуйтесь следующими функциями:


In [93]:
def getCoefficients(model):
    """
    Возвращает отсортированные по абсолютному значению коэффициенты модели
    """
    
    coefs = pd.DataFrame(model.coef_, X_train.columns)
    coefs.columns = ["coef"]
    coefs["abs"] = coefs.coef.apply(np.abs)
    return coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)    
    

def plotCoefficients(model):
    """
    Отрисовывает отсортированные по абсолютному значению коэффициенты модели
    """
    coefs = getCoefficients(model)
    
    plt.figure(figsize=(20, 7))
    coefs.coef.plot(kind='bar')
    plt.grid(True, axis='y')
    plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed');

In [94]:
plotModelResults(lin_reg, X_train, X_test, plot_intervals=True)
plotCoefficients(lin_reg)


Замечательно, интервалы построены, прогнозы достаточно точные, и всё бы хорошо, если бы не одно но. В модели сейчас находится множество признаков и, возможно, какие-то из них лишние, а если их убрать, качество модели серьезно не пострадает. Чтобы убедиться, что у нас есть лишние признаки, постройте тепловую карту корреляций по X_train, поспользовавшись функцией heatmap из библиотеки seaborn:


In [95]:
plt.figure(figsize=(12, 10))
sns.heatmap(X_train);


Действительно, налицо множество корреляций между нашими признаками, и даже видна "сезонность" в этих корреляциях на каждом 24-м лаге. Попробуем регуляризовать нашу модель и убрать несколько лишних переменных.

Обучите на полученных отмасштабированных данных Лассо-регрессию на кросс-валидации (LassoCV), снова передав в качестве параметра cv созданный ранее генератор разбиений.

Постройте график прогнозов с интервалами и убедитесь, что ошибка на тестовой выборке практически не изменилась.


In [96]:
lasso = LassoCV(cv=tscv)
lasso.fit(X_train_scaled, y_train);

In [97]:
plotModelResults(lasso, X_train_scaled, X_test_scaled, plot_intervals=True)


Замечательно, качество нас всё ещё устраивает, интервалы по-прежднему адекватные, посмотрим, сильно ли мы упростили модель. Воспользуйтесь предложенной функцией getCoefficients и скажите, сколько признаков вошли в модель с зануленными весами.


In [98]:
getCoefficients(lasso)


Out[98]:
coef
lag_12 5341.456730
lag_24 3358.843562
lag_47 2852.765783
weekday -2335.154419
lag_14 -1776.472795
lag_25 1335.253935
is_weekend 1271.479455
lag_19 1208.831198
hour 1071.508211
lag_36 -1014.256931
lag_40 -1003.049106
lag_29 926.410445
lag_35 -805.311443
lag_23 799.933271
lag_26 775.350695
lag_18 763.060373
lag_34 -681.241436
lag_37 -586.891345
lag_38 -584.044188
lag_27 573.436411
lag_46 560.249166
lag_33 -532.108098
lag_28 517.795732
lag_13 -450.518763
lag_22 435.256022
lag_42 -396.118203
lag_30 389.584059
lag_41 -329.751723
lag_44 259.617533
lag_43 220.450142
lag_39 -201.267195
lag_17 159.833262
lag_21 17.499581
lag_20 0.000000
lag_45 0.000000
lag_32 -0.000000
lag_16 0.000000
lag_15 -0.000000
lag_31 0.000000

Вопрос 2. Сколько коэффициентов линейной модели равны нулю (с точностью до 6-го знака)?

  • 1
  • 7
  • 5
  • 6

Итак, признаковое пространство хоть немного, но удалось сократить. Что если пойти дальше и сжать наши линейно зависимые признаки в более компактное представление? Для этого воспользуемся методом главных компонент.


In [99]:
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline

def plotPCA():
    """
    График накопленного процента объясненной дисперсии по компонентам
    """
    features = range(pca.n_components_)
    variance = np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)
    plt.figure(figsize=(20, 10))
    plt.bar(features, variance)
    
    # дополнительно отметим уровень, при котором объяснены 95% дисперсии
    plt.hlines(y = 95, xmin=0, xmax=len(features), linestyles='dashed', colors='red')
    
    plt.xlabel('PCA components')
    plt.ylabel('variance')
    plt.xticks(features)
    plt.show()

In [107]:
pca = PCA(0.95)

pca.fit(X_train_scaled)
variances = np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)
plt.plot(variances)
variances


Out[107]:
array([47.91, 68.29, 86.74, 92.15, 94.46, 95.68])

In [64]:
pca.n_components_


Out[64]:
6

Вопрос 3. Какое минимальное число компонент необходимо взять, чтобы объяснить минимум 95% всей дисперсии тренировочного датасета?

  • 4
  • 5
  • 6
  • 7

Снова создайте объект pca, указав при этом оптимальное число компонент (объясняющее минимум 95% изменчивости). Затем создайте две новых переменных – pca_features_train и pca_features_test, присвоив им, соответственно, преобразованные при помощи PCA отмасштабированные датафреймы.


In [80]:
pca = PCA(n_components=6)
pca.fit(X_train_scaled)
pca_features_train = pca.transform(X_train_scaled)
pca_features_test = pca.transform(X_test_scaled)

In [81]:
pca_features_train.shape, pca_features_test.shape


Out[81]:
((1460, 6), (579, 6))

На полученных данных снова обучите простую линейную регрессию и постройте график прогноза.


In [82]:
lin_reg = LinearRegression()
lin_reg.fit(pca_features_train, y_train)
plotModelResults(lin_reg, pca_features_train, pca_features_test, plot_intervals=True)


Вопрос 4. Какая средняя абсолютная ошибка получилась у линейной регрессии, обученной на нескольких главных компонентах?

  • 5140
  • 4917
  • 6719
  • 5434